stanEx1
execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}ex1 normal distribution
ex1.stan
data{
int N;
vector[N] y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~normal(m,s);
}mdl=cmdstan_model('./ex1.stan')
y=rnorm(10,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.662 1.042 2.136 1.901 2.546 2.904par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y') 
data=list(N=length(y),y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -1357.56
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -3.2695 4.48716e-05 1.26718e-06 1 1 23
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.3 seconds.mle## variable estimate
## lp__ -3.27
## m 1.90
## s 0.84mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -4.51 -4.17 1.15 0.78 -6.75 -3.46 1.00 802 856
## m 1.93 1.93 0.34 0.30 1.39 2.48 1.01 890 693
## s 1.04 0.98 0.32 0.24 0.68 1.60 1.00 1194 1099mcmc$metadata()$stan_variables## [1] "lp__" "m" "s"mcmc$metadata()$model_params## [1] "lp__" "m" "s"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -4.51 -4.17 1.15 0.78 -6.75 -3.46 1.00 802 856
## m 1.93 1.93 0.34 0.30 1.39 2.48 1.01 890 693
## s 1.04 0.98 0.32 0.24 0.68 1.60 1.00 1194 1099
ex2 poisson distribution
ex2
y=rpois(20,3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 1.00 2.50 2.65 4.00 6.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=length(y),y=y)ex2-1.stan
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> m;
real<lower=0> s;
}
model{
y~normal(m,s);
}
generated quantities{
vector[N] y1;
for(i in 1:N)
y1[i]=normal_rng(m,s);
}# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -38.2607
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 7 -20.3939 0.000406511 0.000836652 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -20.39
## m 2.65
## s 1.68
## y1[1] 3.39
## y1[2] 0.92
## y1[3] -0.01
## y1[4] 3.33
## y1[5] 1.89
## y1[6] 2.91
## y1[7] 3.97
## y1[8] 2.69
## y1[9] 2.26
## y1[10] 2.76
## y1[11] 4.05
## y1[12] 3.32
## y1[13] 5.38
## y1[14] 3.00
## y1[15] 3.31
## y1[16] 6.28
## y1[17] 2.42
## y1[18] 4.76
## y1[19] 1.93
## y1[20] -1.42mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -19.95 -19.60 1.04 0.76 -22.10 -18.94 1.00 757 1110
## m 2.66 2.66 0.41 0.41 1.99 3.33 1.01 1420 1128
## s 1.85 1.80 0.33 0.31 1.39 2.49 1.00 1180 952
## y1[1] 2.66 2.64 1.95 1.81 -0.58 5.83 1.00 1658 1904
## y1[2] 2.63 2.70 1.97 1.84 -0.62 5.85 1.00 1771 1317
## y1[3] 2.66 2.65 1.94 1.90 -0.49 5.84 1.00 2050 1747
## y1[4] 2.67 2.63 1.88 1.82 -0.36 5.72 1.00 2008 2059
## y1[5] 2.62 2.59 1.93 1.88 -0.48 5.82 1.00 1916 1822
## y1[6] 2.68 2.66 1.94 1.85 -0.50 5.92 1.00 1980 1999
## y1[7] 2.66 2.71 1.88 1.76 -0.40 5.82 1.00 2004 1933
## y1[8] 2.69 2.61 1.90 1.80 -0.28 5.94 1.00 1867 2007
## y1[9] 2.63 2.63 1.94 1.83 -0.54 5.83 1.00 1984 1969
## y1[10] 2.70 2.73 1.89 1.74 -0.51 5.77 1.00 1995 1876
## y1[11] 2.71 2.70 1.87 1.84 -0.35 5.82 1.00 2019 1854
## y1[12] 2.62 2.67 1.86 1.87 -0.44 5.61 1.00 2097 1854
## y1[13] 2.60 2.60 1.92 1.82 -0.53 5.86 1.00 1924 1901
## y1[14] 2.68 2.66 1.89 1.81 -0.42 5.73 1.00 2003 1938
## y1[15] 2.71 2.75 1.98 1.96 -0.53 6.10 1.00 2085 1972
## y1[16] 2.59 2.62 1.95 1.89 -0.56 5.79 1.00 1896 1974
## y1[17] 2.65 2.58 1.91 1.85 -0.51 5.83 1.00 1958 2052
## y1[18] 2.66 2.68 1.94 1.90 -0.47 5.71 1.00 2017 1919
## y1[19] 2.65 2.60 1.94 1.87 -0.58 5.89 1.00 1902 1954
## y1[20] 2.66 2.68 1.97 1.84 -0.69 5.91 1.00 1950 1971
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
ex2-2.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0> l;
}
model{
y~poisson(l);
}
generated quantities{
array[N] int y1;
for(i in 1:N)
y1[i]=poisson_rng(l);
}mdl=cmdstan_model('./ex2-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -1.52152
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 -1.34834 1.36669e-05 7.81676e-06 0.99 0.99 6
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -1.35
## l 2.65
## y1[1] 4.00
## y1[2] 2.00
## y1[3] 3.00
## y1[4] 2.00
## y1[5] 2.00
## y1[6] 0.00
## y1[7] 2.00
## y1[8] 2.00
## y1[9] 3.00
## y1[10] 1.00
## y1[11] 2.00
## y1[12] 2.00
## y1[13] 5.00
## y1[14] 2.00
## y1[15] 0.00
## y1[16] 2.00
## y1[17] 5.00
## y1[18] 5.00
## y1[19] 1.00
## y1[20] 3.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -0.84 -0.57 0.68 0.28 -2.16 -0.37 1.00 1007 1383
## l 2.70 2.68 0.36 0.35 2.14 3.32 1.00 812 1179
## y1[1] 2.70 3.00 1.66 1.48 0.00 6.00 1.00 1813 1763
## y1[2] 2.67 2.00 1.69 1.48 0.00 6.00 1.00 1971 1892
## y1[3] 2.73 3.00 1.69 1.48 0.00 6.00 1.00 1894 1902
## y1[4] 2.68 2.00 1.64 1.48 0.00 6.00 1.00 1890 1918
## y1[5] 2.67 2.00 1.69 1.48 0.00 6.00 1.00 1933 2122
## y1[6] 2.72 3.00 1.69 1.48 0.00 6.00 1.00 1744 1890
## y1[7] 2.66 3.00 1.61 1.48 0.00 5.00 1.00 2048 2033
## y1[8] 2.70 3.00 1.70 1.48 0.00 6.00 1.00 2075 1933
## y1[9] 2.73 3.00 1.71 1.48 0.00 6.00 1.00 1883 1978
## y1[10] 2.67 2.00 1.65 1.48 0.00 6.00 1.00 1697 1875
## y1[11] 2.67 3.00 1.65 1.48 0.00 6.00 1.00 1952 1916
## y1[12] 2.66 2.00 1.67 1.48 0.00 6.00 1.00 1851 1833
## y1[13] 2.68 3.00 1.70 1.48 0.00 6.00 1.00 1886 1936
## y1[14] 2.74 3.00 1.66 1.48 0.00 6.00 1.00 1888 1829
## y1[15] 2.70 2.00 1.73 1.48 0.00 6.00 1.00 1923 1924
## y1[16] 2.64 2.00 1.66 1.48 0.00 6.00 1.00 1797 1824
## y1[17] 2.75 3.00 1.68 1.48 0.00 6.00 1.00 1710 1770
## y1[18] 2.66 2.00 1.66 1.48 0.00 6.00 1.00 2051 1845
## y1[19] 2.70 2.00 1.70 1.48 0.00 6.00 1.00 1799 2001
## y1[20] 2.65 2.00 1.64 1.48 0.00 6.00 1.00 1986 2102
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
ex3
difference test
ex3.stan
data{
int N;
vector[N] a;
vector[N] b;
}
parameters{
real ma;
real<lower=0> sa;
real mb;
real<lower=0> sb;
}
model{
a~normal(ma,sa);
b~normal(mb,sb);
}
generated quantities{
real d;
d=mb-ma;
}n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)
mdl=cmdstan_model('./ex3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -6122.02
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 29 -35.2022 0.000122699 0.000172573 1 1 44
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -35.20
## ma 10.16
## sa 1.10
## mb 11.21
## sb 1.94
## d 1.05mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -36.45 -36.15 1.42 1.27 -39.16 -34.79 1.00 792 1157
## ma 10.17 10.16 0.27 0.25 9.73 10.61 1.00 1388 1303
## sa 1.21 1.18 0.21 0.20 0.92 1.61 1.00 2003 1484
## mb 11.21 11.22 0.47 0.46 10.41 11.96 1.01 822 833
## sb 2.14 2.09 0.38 0.35 1.61 2.84 1.00 2235 1408
## d 1.04 1.04 0.54 0.53 0.14 1.92 1.00 955 887
d=mcmc$draws('d')
d## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
##
## chain
## iteration 1 2 3 4
## 1 0.94 0.50 0.96 1.51
## 2 1.27 0.69 0.77 1.43
## 3 1.44 0.63 0.39 0.45
## 4 1.44 0.55 0.64 2.80
## 5 1.54 1.75 0.76 2.04
##
## # ... with 495 more iterationsmean(d>0)## [1] 0.971ex4
ex4-1.stan
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)
data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -291608
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 36 -50.5368 0.00179027 0.000612657 0.9903 0.9903 84
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -50.54
## b0 8.72
## b1 3.00
## s 7.59mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -50.03 -49.72 1.21 1.02 -52.31 -48.68 1.00 510 474
## b0 8.68 8.58 4.16 3.92 1.96 15.52 1.01 249 260
## b1 3.00 3.00 0.07 0.06 2.89 3.11 1.01 305 375
## s 8.63 8.47 1.57 1.49 6.45 11.37 1.00 931 1045
b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)
ex4-2.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N1] m;
vector[N1] y1;
for(i in 1:N1){
m[i]=b0+b1*x1[i];
y1[i]=normal_rng(m[i],s);
}
}# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex4-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -3431.79
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 30 -50.5368 0.0262534 0.000427242 0.9908 0.9908 54
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -50.54
## b0 8.72
## b1 3.00
## s 7.59
## m[1] 48.61
## m[2] 76.03
## m[3] 103.44
## m[4] 130.86
## m[5] 158.27
## m[6] 185.69
## m[7] 213.10
## m[8] 240.52
## m[9] 267.93
## m[10] 295.35
## y1[1] 44.85
## y1[2] 80.67
## y1[3] 102.04
## y1[4] 131.69
## y1[5] 150.33
## y1[6] 190.10
## y1[7] 215.18
## y1[8] 232.94
## y1[9] 272.21
## y1[10] 301.88mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -50.09 -49.75 1.30 1.03 -52.55 -48.71 1.01 421 603
## b0 9.23 9.18 4.42 4.31 1.95 16.95 1.01 216 258
## b1 2.99 2.99 0.07 0.07 2.87 3.10 1.01 256 379
## s 8.61 8.39 1.50 1.42 6.57 11.28 1.00 885 798
## m[1] 49.03 48.89 3.61 3.48 43.02 55.37 1.01 219 265
## m[2] 76.38 76.24 3.09 2.95 71.24 81.70 1.01 227 274
## m[3] 103.73 103.63 2.62 2.58 99.56 108.13 1.01 248 302
## m[4] 131.09 131.06 2.25 2.22 127.62 134.66 1.01 303 309
## m[5] 158.44 158.44 2.01 1.96 155.24 161.59 1.00 476 547
## m[6] 185.79 185.82 1.96 1.86 182.63 188.85 1.00 1154 1260
## m[7] 213.15 213.12 2.11 2.01 209.77 216.52 1.00 2077 1397
## m[8] 240.50 240.54 2.43 2.36 236.65 244.47 1.00 1758 1322
## m[9] 267.85 267.83 2.86 2.77 263.29 272.65 1.00 1006 1329
## m[10] 295.21 295.17 3.36 3.23 289.88 300.81 1.00 699 1185
## y1[1] 49.35 49.14 9.48 9.42 34.30 64.92 1.00 1065 1494
## y1[2] 76.71 76.80 9.06 8.75 61.65 91.33 1.00 1288 1466
## y1[3] 103.47 103.65 9.11 8.83 88.77 118.28 1.00 1355 1780
## y1[4] 130.97 131.04 8.78 8.47 116.64 145.04 1.00 1730 1730
## y1[5] 159.18 159.02 9.04 8.59 143.77 174.21 1.00 1909 1889
## y1[6] 185.55 185.52 8.91 8.76 171.16 200.23 1.00 2092 1887
## y1[7] 213.14 213.06 8.90 8.75 198.91 227.82 1.00 1653 1816
## y1[8] 240.40 240.56 9.16 8.88 225.17 255.09 1.00 1820 1894
## y1[9] 268.00 267.95 9.33 9.00 252.95 283.27 1.00 1960 1932
## y1[10] 295.37 295.19 9.37 8.58 280.13 310.57 1.00 1565 1649
m=mcmc$draws('m')
summary(m)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 m[1] 49.0 48.9 3.61 3.48 43.0 55.4 1.01 220. 265.
## 2 m[2] 76.4 76.2 3.09 2.95 71.2 81.7 1.01 228. 274.
## 3 m[3] 104. 104. 2.62 2.58 99.6 108. 1.01 248. 302.
## 4 m[4] 131. 131. 2.25 2.22 128. 135. 1.01 304. 310.
## 5 m[5] 158. 158. 2.01 1.96 155. 162. 1.00 476. 548.
## 6 m[6] 186. 186. 1.96 1.86 183. 189. 1.00 1155. 1260.
## 7 m[7] 213. 213. 2.11 2.01 210. 217. 1.00 2078. 1398.
## 8 m[8] 240. 241. 2.43 2.36 237. 244. 1.00 1758. 1322.
## 9 m[9] 268. 268. 2.86 2.77 263. 273. 1.00 1006. 1330.
## 10 m[10] 295. 295. 3.36 3.23 290. 301. 1.00 700. 1186.smm=summary(m)
y1=mcmc$draws('y1')
summary(y1)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 y1[1] 49.3 49.1 9.48 9.42 34.3 64.9 1.00 1065. 1494.
## 2 y1[2] 76.7 76.8 9.06 8.75 61.7 91.3 1.00 1289. 1466.
## 3 y1[3] 103. 104. 9.11 8.83 88.8 118. 1.00 1356. 1780.
## 4 y1[4] 131. 131. 8.78 8.47 117. 145. 1.00 1731. 1730.
## 5 y1[5] 159. 159. 9.04 8.59 144. 174. 1.00 1909. 1890.
## 6 y1[6] 186. 186. 8.91 8.76 171. 200. 1.00 2093. 1888.
## 7 y1[7] 213. 213. 8.90 8.75 199. 228. 1.00 1653. 1817.
## 8 y1[8] 240. 241. 9.16 8.88 225. 255. 1.00 1821. 1894.
## 9 y1[9] 268. 268. 9.33 9.00 253. 283. 1.00 1961. 1932.
## 10 y1[10] 295. 295. 9.37 8.58 280. 311. 1.00 1565. 1650.smy=summary(y1)
xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)
par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=m),xy,col='black')+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')
ex4-3.stan
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m;
vector[N] y1;
for(i in 1:N){
m[i]=b0+b1*x[i];
y1[i]=normal_rng(m[i],s);
}
}data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-3.stan')
mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -50.05 -49.65 1.36 0.99 -52.86 -48.66 1.02 346 218
## b0 9.02 8.79 4.34 3.79 2.47 16.12 1.05 108 105
## b1 2.99 2.99 0.07 0.06 2.88 3.10 1.04 153 147
## s 8.66 8.43 1.62 1.40 6.42 11.64 1.01 667 659
## m[1] 69.42 69.25 3.14 2.85 64.69 74.50 1.05 102 88
## m[2] 64.71 64.54 3.23 2.92 59.88 69.98 1.05 103 85
## m[3] 131.00 130.86 2.16 2.01 127.67 134.50 1.03 150 184
## m[4] 48.85 48.68 3.54 3.19 43.53 54.59 1.05 106 106
## m[5] 245.84 245.81 2.38 2.28 241.89 249.83 1.00 1551 1377
## m[6] 257.86 257.83 2.57 2.50 253.58 262.07 1.00 1186 1137
## m[7] 293.54 293.53 3.19 3.13 288.23 298.56 1.01 669 826
## m[8] 144.56 144.42 2.02 1.85 141.43 147.80 1.03 197 260
## m[9] 256.76 256.73 2.55 2.48 252.50 260.94 1.00 1214 1159
## m[10] 246.59 246.57 2.40 2.30 242.62 250.59 1.00 1526 1298
## m[11] 169.93 169.89 1.86 1.71 166.96 173.04 1.01 472 482
## m[12] 295.20 295.18 3.22 3.16 289.83 300.26 1.01 655 826
## m[13] 292.87 292.86 3.18 3.12 287.59 297.89 1.01 675 801
## m[14] 192.33 192.35 1.87 1.72 189.25 195.43 1.00 1104 1300
## m[15] 112.33 112.19 2.42 2.25 108.73 116.20 1.04 111 135
## m[16] 272.50 272.47 2.81 2.74 267.74 277.05 1.01 894 978
## m[17] 128.08 127.95 2.20 2.08 124.71 131.69 1.03 144 178
## m[18] 94.95 94.80 2.69 2.50 90.94 99.19 1.04 101 91
## m[19] 91.03 90.90 2.76 2.54 86.90 95.34 1.04 100 91
## m[20] 115.46 115.33 2.37 2.21 111.90 119.22 1.04 115 135
## y1[1] 69.53 69.61 9.18 8.84 54.45 84.47 1.01 1087 1936
## y1[2] 65.05 65.05 9.19 9.05 50.90 80.39 1.01 1024 974
## y1[3] 131.00 130.90 8.85 8.67 116.54 145.07 1.00 1691 1480
## y1[4] 48.90 48.73 9.48 8.70 33.72 64.47 1.01 568 1451
## y1[5] 245.62 245.49 9.17 8.83 230.45 260.66 1.00 2007 1564
## y1[6] 258.37 258.37 9.11 8.48 243.19 273.56 1.00 1879 2053
## y1[7] 293.45 293.51 9.40 9.01 278.12 308.61 1.00 1760 1566
## y1[8] 144.66 144.76 9.00 8.86 129.81 159.51 1.00 1630 1916
## y1[9] 256.84 257.19 9.18 8.79 240.99 271.41 1.00 1988 1916
## y1[10] 246.74 246.88 9.13 8.69 231.05 261.44 1.00 1897 1933
## y1[11] 169.82 169.56 9.22 8.55 154.84 185.25 1.00 2001 1618
## y1[12] 295.01 295.26 9.21 8.71 279.94 310.23 1.00 1875 1841
## y1[13] 293.37 293.38 9.13 8.66 278.28 307.77 1.00 1610 1916
## y1[14] 192.45 192.46 8.90 8.68 178.07 206.77 1.00 2046 1810
## y1[15] 112.51 112.53 9.14 8.72 98.03 127.82 1.00 1517 1890
## y1[16] 272.52 272.34 9.28 8.46 257.61 287.76 1.00 1646 1816
## y1[17] 128.04 128.01 9.12 8.68 113.11 142.53 1.00 1619 1599
## y1[18] 95.08 95.28 9.06 8.86 80.06 109.78 1.00 1412 1673
## y1[19] 90.73 90.72 9.20 9.14 75.93 106.04 1.00 1149 1785
## y1[20] 115.27 115.36 9.18 8.63 99.96 129.84 1.00 1607 1736
y1=mcmc$draws('y1')
smy=summary(y1)
par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
estimate correlation ex4-41.stan
data{
int N;
array[N] vector[2] y;
}
parameters{
vector[2] m;
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
y~multi_normal(m,cv);
}ex4-42.stan
data{
int N;
matrix[2,2] dp;
}
parameters{
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
dp~wishart(N-1,cv);
}ex4-43.stan
data{
int N;
int M;
matrix[M,M] dp;
}
parameters{
vector<lower=0>[M] s;
corr_matrix[M] r;
}
transformed parameters{
cov_matrix[M] cv;
cv=quad_form_diag(r,s);
}
model{
dp~wishart(N-1,cv);
}ex4-44.stan
data{
int N;
int K;
array[N] vector[K] y;
}
parameters{
vector[K] m;
cov_matrix[K] cov;
}
model{
y~multi_normal(m,cov);
}library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)## [,1] [,2]
## [1,] 0.867 0.538
## [2,] 0.538 0.784cor(y)## [,1] [,2]
## [1,] 1.000 0.653
## [2,] 0.653 1.000plot(y)
#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')
data=list(N=n,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -35.5849
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -9.54568 2.07811e-05 0.000529216 0.7427 0.7427 18
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -9.55
## m[1] 0.07
## m[2] 0.25
## s[1] 0.91
## s[2] 0.86
## r 0.65
## cv[1,1] 0.82
## cv[2,1] 0.51
## cv[1,2] 0.51
## cv[2,2] 0.74mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -13.71 -13.37 1.72 1.48 -16.90 -11.59 1.01 825 1159
## m[1] 0.07 0.07 0.24 0.23 -0.32 0.45 1.00 1295 1239
## m[2] 0.24 0.25 0.24 0.22 -0.14 0.62 1.00 1355 1330
## s[1] 1.01 0.98 0.18 0.16 0.76 1.34 1.00 1284 1300
## s[2] 0.96 0.93 0.18 0.16 0.71 1.30 1.00 1183 1251
## r 0.61 0.63 0.15 0.14 0.32 0.81 1.02 443 679
## cv[1,1] 1.05 0.97 0.40 0.31 0.58 1.79 1.00 1284 1300
## cv[2,1] 0.62 0.55 0.31 0.26 0.24 1.21 1.01 550 644
## cv[1,2] 0.62 0.55 0.31 0.26 0.24 1.21 1.01 550 644
## cv[2,2] 0.96 0.87 0.39 0.29 0.51 1.69 1.00 1183 1251
#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n
mdl=cmdstan_model('./ex4-42.stan')
data=list(N=n,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -324.231
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 14 -10.043 9.20225e-05 0.000306709 1 1 18
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -10.04
## s[1] 0.93
## s[2] 0.89
## r 0.65
## cv[1,1] 0.87
## cv[2,1] 0.54
## cv[1,2] 0.54
## cv[2,2] 0.78mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -13.04 -12.72 1.32 1.05 -15.60 -11.63 1.01 830 687
## s[1] 1.01 0.99 0.19 0.17 0.76 1.34 1.00 1395 1240
## s[2] 0.95 0.93 0.17 0.16 0.72 1.26 1.00 1381 1160
## r 0.60 0.62 0.16 0.14 0.30 0.80 1.01 455 506
## cv[1,1] 1.06 0.98 0.44 0.33 0.58 1.80 1.00 1395 1240
## cv[2,1] 0.60 0.55 0.31 0.25 0.23 1.14 1.01 564 603
## cv[1,2] 0.60 0.55 0.31 0.25 0.23 1.14 1.01 564 603
## cv[2,2] 0.94 0.87 0.36 0.30 0.52 1.59 1.00 1381 1160
#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')
m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -1230.72
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -10.043 5.06124e-05 1.12298e-05 1 1 21
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -10.04
## s[1] 0.93
## s[2] 0.89
## r[1,1] 1.00
## r[2,1] 0.65
## r[1,2] 0.65
## r[2,2] 1.00
## cv[1,1] 0.87
## cv[2,1] 0.54
## cv[1,2] 0.54
## cv[2,2] 0.78mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.30 -11.97 1.28 1.04 -14.79 -10.92 1.00 812 1252
## s[1] 1.01 0.98 0.18 0.17 0.76 1.34 1.00 953 1146
## s[2] 0.96 0.94 0.17 0.16 0.73 1.27 1.00 1121 1241
## r[1,1] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## r[2,1] 0.61 0.63 0.15 0.14 0.33 0.81 1.00 735 1032
## r[1,2] 0.61 0.63 0.15 0.14 0.33 0.81 1.00 735 1032
## r[2,2] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## cv[1,1] 1.05 0.96 0.41 0.33 0.58 1.80 1.00 953 1146
## cv[2,1] 0.62 0.56 0.30 0.25 0.25 1.20 1.00 684 747
## cv[1,2] 0.62 0.56 0.30 0.25 0.25 1.20 1.00 684 747
## cv[2,2] 0.95 0.88 0.37 0.30 0.53 1.63 1.00 1121 1241
mdl=cmdstan_model('./ex4-44.stan')
k=2
data=list(N=n,K=k,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -1110.98
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -9.54568 1.35409e-05 0.000670983 0.5952 0.5952 24
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -9.55
## m[1] 0.07
## m[2] 0.25
## cov[1,1] 0.82
## cov[2,1] 0.51
## cov[1,2] 0.51
## cov[2,2] 0.74mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.04 -11.70 1.77 1.67 -15.39 -9.82 1.01 827 1236
## m[1] 0.08 0.07 0.26 0.24 -0.35 0.51 1.00 1062 1160
## m[2] 0.25 0.25 0.24 0.22 -0.14 0.64 1.00 1095 1033
## cov[1,1] 1.27 1.14 0.55 0.42 0.65 2.35 1.00 1330 1338
## cov[2,1] 0.78 0.69 0.44 0.35 0.26 1.61 1.00 926 891
## cov[1,2] 0.78 0.69 0.44 0.35 0.26 1.61 1.00 926 891
## cov[2,2] 1.14 1.03 0.50 0.39 0.56 2.07 1.00 1209 1127
